lapop12 = readRDS('data/lapop12.rds')

#### Imputing #### 

lapop12$soc.cl_miss <- as.numeric(lapop12$soc.cl_bin >= 0, 1,0)
lapop12$soc.cl_miss[is.na(lapop12$soc.cl_miss)] <- 0 

lapop12$torture_miss <- as.numeric(lapop12$torture_bin >= 0, 1,0)
lapop12$torture_miss[is.na(lapop12$torture_miss)] <- 0 

lapop12$trial_miss <- as.numeric(lapop12$fair_trial > 0, 1,0)
lapop12$trial_miss[is.na(lapop12$trial_miss)] <- 0 

lapop_df2 <- lapop12 %>% 
  group_by(ur, estratopri, estratosec) %>%
  summarize(soc.cl_mean  = mean(soc.cl_bin, na.rm = T), 
            torture_mean = mean(torture_bin, na.rm = T), 
            fair_trial_mean = mean(fair_trial, na.rm = T))

lapop_df <- merge(lapop12, lapop_df2, by = c("ur", "estratopri", "estratosec"))
lapop_df$soc.cl_impute <- ifelse(lapop_df$soc.cl_miss==1, lapop_df$soc.cl_bin, lapop_df$soc.cl_mean)
lapop_df$torture_impute <- ifelse(lapop_df$torture_miss==1, lapop_df$torture_bin, lapop_df$torture_mean)
lapop_df$fair_trial_impute <- ifelse(lapop_df$trial_miss==1, lapop_df$fair_trial, lapop_df$fair_trial_mean)


m1 <- lm_robust(torture_impute ~ trial + torture_miss, data = lapop_df, 
                fixed_effects = estratopri + tamano + ur, 
                se_type = 'stata')

m2 <- lm_robust(soc.cl_impute ~ trial + soc.cl_miss, data = lapop_df, 
                fixed_effects = estratopri + tamano + ur, 
                se_type = 'stata')

m3 <- lm_robust(fair_trial_impute ~ trial + trial_miss, data = lapop_df, 
                fixed_effects = estratopri + tamano + ur, 
                se_type = 'stata')

coef.name <- c('Torture', 'Soc. Cleanse', 'Fair Trial')
y <- c(summary(m1)$coefficients[1,1], summary(m2)$coefficients[1,1], summary(m3)$coefficients[1,1])  
U <- c(summary(m1)$coefficients[1,6], summary(m2)$coefficients[1,6], summary(m3)$coefficients[1,6])  
L <- c(summary(m1)$coefficients[1,5], summary(m2)$coefficients[1,5], summary(m3)$coefficients[1,5])  
pv <- c(m1$p.value[1], m2$p.value[1],  
        m3$p.value[1])
plot.df <- data.frame(coef.name, y, U, L, pv)
plot.df$y2 <- round(plot.df$y, 2)
plot.df$y2 <- ifelse(plot.df$pv < .05 & plot.df$pv > .01, paste0(plot.df$y2, "*"), plot.df$y2)
plot.df$y2 <- ifelse(plot.df$pv < .01 & plot.df$pv > .001, paste0(plot.df$y2, "**"), plot.df$y2)
plot.df$y2 <- ifelse(plot.df$pv < .001, paste0(plot.df$y2, "***"), plot.df$y2)


impute <- ggplot(plot.df, aes(x = coef.name, y = y)) +
  geom_point(size = 1) +
  geom_errorbar(aes(ymax = U, ymin = L), width=.1,
                position=position_dodge(.9)) + 
  geom_hline(yintercept = 0) + 
  ggtitle("Imputed Block Mean") + 
  xlab("Outcome") +
  ylab("Partial Derivative") +
  annotate("label", x = coef.name, y = y, label = plot.df$y2, 
           hjust=.5, vjust=1, family="Times") +
  theme_tufte()

impute

ggsave("fig-out/impute.pdf", width = 8)

